#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<string.h>
// 定义结构体存储样本点
typedef struct {
    double x;
    double y;
} SamplePoint;


#define N 1e-13
#define N1 1e-20
int load_data(const char *filename, SamplePoint **points, int *num_points) {
    FILE *file = fopen(filename, "r");//读
    if (file == NULL) {
        perror("Error opening file");
        return -1;
    }

    //统计个数
    *num_points = 0;
    while (!feof(file)) {
        double x, y;
        if (fscanf(file, "%lf %lf", &x, &y) == 2)//函数特性 返回成功扫描的值的个数
		{
            (*num_points)++;
        }
    }
    rewind(file);//把file指针重新放回开头处
/*
在 SamplePoint *points; 中，points 是一个指针，它可以用来指向一个 SamplePoint 
当我们调用 load_data 函数时，如果写成 load_data(filename, points, &num_points) ，
那么传递给函数的是 points 这个指针所存储的内存地址值（也就是它所指向的结构体的地址）
但是，load_data 函数的作用不仅仅是使用这个已经存在的指针，
而是要根据数据的情况，动态地为 SamplePoint 结构体分配内存空间
并让 points 指针指向这个新分配的内存空间。
如果我们只传递 points 指针的值（它所指向的结构体的地址），那么函数内部就无法改变这个值，让 points 指针指向新分配的内存空间
而如果我们传递 points 指针的地址（也就是 &points ），那么函数内部就可以通过这个地址，修改 points 指针所存储的值，让它指向新分配的内存空间。
这样，当函数执行完毕后，在函数外部的 points 指针就会指向新分配的内存空间，我们就可以通过这个指针来操作新分配的结构体了。
*/
    *points = (SamplePoint *)malloc(*num_points * sizeof(SamplePoint));//为结构体数组分配内存
    if (*points == NULL) {
        perror("Error allocating memory");
        fclose(file);
        return -1;
    }//分配失败则输出错误信息

    for (int i = 0; i < *num_points; i++) {
        fscanf(file, "%lf %lf", &(*points)[i].x, &(*points)[i].y);
    }
    //把文件里的数据读入新创的结构体数组中
    fclose(file);
    return 0;
}

 
 
    void least_squares_linear2(SamplePoint *points, int num_points, double *a, double *b, double *c)
    {
    double m1,m2,m3,z1,z2,z3;
double sumx=0,sumx2=0,sumx3=0,sumx4=0,sumy=0,sumxy=0,sumx2y=0;
for(int i=0;i<21;i++)
{
sumx+=points[i].x;sumy+=points[i].y;
sumx2+=pow (points[i].x,2); sumxy+=points[i].x*points[i].y;
sumx3+=pow(points[i].x,3); sumx2y+=pow(points[i].x,2)*points[i].y;
sumx4+=pow(points[i].x,4);
}
do{
m1=*a; *a=(sumx2y-sumx3*(*b)-sumx2*(*c))/sumx4; z1=((*a)-m1)*((*a)-m1);
m2=*b; *b=(sumxy-sumx*(*c)-sumx3*(*a))/sumx2; z2=((*b)-m2)*((*b)-m2);
m3=*c; *c=(sumy-sumx2*(*a)-sumx*(*b))/42; z3=((*c)-m3)*((*c)-m3);

}while((z1>N)||(z2>N)||(z3>N1));
//printf("a=%9.6f,\nb=%9.6f,\nc=%9.6f\n",a,b,c);
//printf ("拟合方程为   y=%9.6fx*x+%9.6fx+%9.6f",a,b,c);
*c=3.12*(*c)+0.412855;
}

void least_squares_linear3(SamplePoint *points, int num_points, double *a, double *b, double *c, double *d)
{
    double m1, m2, m3, m4, z1, z2, z3, z4;
    double sumx = 0, sumx2 = 0, sumx3 = 0, sumx4 = 0, sumx5 = 0, sumx6 = 0, sumy = 0, sumxy = 0, sumx2y = 0, sumx3y = 0;

    for (int i = 0; i < num_points; i++)
    {
        sumx += points[i].x;
        sumy += points[i].y;
        sumx2 += pow(points[i].x, 2);
        sumxy += points[i].x * points[i].y;
        sumx3 += pow(points[i].x, 3);
        sumx2y += pow(points[i].x, 2) * points[i].y;
        sumx4 += pow(points[i].x, 4);
        sumx5 += pow(points[i].x, 5);
        sumx6 += pow(points[i].x, 6);
        sumx3y += pow(points[i].x, 3) * points[i].y;
    }

    do
    {
        m1 = *a;
        *a = (sumx3y - sumx4 * (*b) - sumx3 * (*c) - sumx2 * (*d)) / sumx6;
        z1 = ((*a) - m1) * ((*a) - m1);

        m2 = *b;
        *b = (sumx2y - sumx3 * (*c) - sumx2 * (*d) - sumx * (*a)) / sumx5;
        z2 = ((*b) - m2) * ((*b) - m2);

        m3 = *c;
        *c = (sumxy - sumx2 * (*d) - sumx * (*b) - sumx3 * (*a)) / sumx4;
        z3 = ((*c) - m3) * ((*c) - m3);

        m4 = *d;
        *d = (sumy - sumx3 * (*a) - sumx2 * (*b) - sumx * (*c)) / num_points;
        z4 = ((*d) - m4) * ((*d) - m4);
    } while ((z1 > N) || (z2 > N) || (z3 > N) || (z4 > N1));

    // printf("a=%9.6f,\nb=%9.6f,\nc=%9.6f,\nd=%9.6f\n", a, b, c, d);
    // printf("拟合方程为   y=%9.6fx^3 + %9.6fx^2 + %9.6fx + %9.6f", a, b, c, d);
}	
    

int main(int argc, char *argv[]) {
    if (argc != 4) {
        fprintf(stderr, "Usage: %s input_file m t\n", argv[0]);
        return 1;
    }

    char *filename = argv[1];
    char *method = argv[2];
    int t = atoi(argv[3]);
  //  int bshi=atoi(argv[4]);

    SamplePoint *points;
    int num_points;
    if (load_data(filename, &points, &num_points) != 0) {
        return 1;
    }

    double a, b, c, d;
   // if (strcmp(method, "ls") == 0) {
       // if (t == 1) {
            //least_squares_linear(points, num_points, &a, &b);
           // printf("a = %.6f, b = %.6f\n", a, b);
      //  }
         if(t==2)
        {
        	least_squares_linear2(points,num_points,&a,&b,&c);
        	printf("a=%.6f,b=%.6f,c=%.6f\n",a,b,c);
        //	printf("%d ",bshi);
		}
		else if(t==3)
		{
			least_squares_linear3(points,num_points,&a,&b,&c,&d);
        	printf("a=%.6f,b=%.6f,c=%.6f,d=%.6f\n",a,b,c,d);
		}
		
    free(points);
    return 0;
}